Talvolta si è soliti rappresentare un fenomeno esaminando il solo impatto delle variabili esplicative prese singolarmente. Questo approccio risulta più semplice e facilmente interpretabile, ma non sempre coerente alla realtà. In molte situazioni è possibile migliorare la capacità previsiva del modello considerando l’effetto congiunto di due o più predittori.
Per poter meglio cogliere l’importanza delle interazioni si consideri il seguente esempio:
“si supponga di voler massimizzare la resa di un campo di grano utilizzando acqua e fertilizzante. Qualora si utilizzasse la quantità ottima di fertilizzante ma senza usare acqua il campo non renderebbe nulla poiché le piante non potrebbero crescere senza acqua. Se, invece, si utilizzasse la quantità ottima d’acqua ma senza fertilizzante la resa del campo sarebbe buona, ma non ottimale, poiché le piante necessiterebbero comunque del fertilizzante.Per massimizzare il rendimento ottenibile dal campo è necessario combinare acqua e fertilizzante nella quantità ottimale poiché questi due fattori interagiscono tra loro.”
Da questo esempio si comprende che, talvolta, per poter rappresentare correttamente un fenomeno è necessario introdurre un’interazione tra le variabili coinvolte.
Formalmente, si può dire che due predittori interagiscono se il loro effetto congiunto sulla variabile dipendente è diverso da quello che ci si aspetta valutandoli singolarmente.
Si consideri un modello lineare, un’interazione è esprimibile come segue:
\[ y= \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + \epsilon \] dove,
\(\beta_0\) l’intercetta, corrisponde al valore medio della variabile dipendente, y, quando \(x_1\) e \(x_2\) sono entrambe a zero.
\(\beta_1\) rappresenta la variazione della variabile dipendente in corrispondenza della variazione unitaria di \(x_1\) se \(x_2\) è nulla.
\(\beta_2\) rappresenta la variazione di \(y\) quando \(x_2\) subisce una variazione unitaria e \(x_1\) è nulla.
\(\beta_3\) è il coefficiente di regressione corrispondente al termine di interazione tra \(x_1\) e \(x_2\) e rappresenta il modo in cui le due variabili influenzano congiuntamente \(y\).
\(\epsilon\) è la componente stocastica del modello, essa rappresenta la variazione di y che non può essere spiegata dalla componente deterministica dell’equazione.
Si consideri il seguente modello:
\[modello\ di\ interesse:\ y_t=\beta_1+\beta_2x_t+\epsilon_t\ ,\ t=1,...,T\] Qualora si fosse interessati a incorporare informazioni qualitative è necessario ricorrere all’utilizzo di una variabile dummy.
\[d_t=\begin{cases} 0 &se\ t\in A \\1&se\ t\notin A \end{cases}\] \[modello\ con\ interazione:\ y_t=\beta_1+\beta_2x_t+\beta_3d_tx_t+\epsilon_t\ ,\ t=1,...,T\] Si noti che \(\beta_3\) non è un effetto marginale, esso è interpretabile come:
\(d_t =0 \rightarrow E(Y_t)=\beta_1+\beta_2x_t\)
\(d_t =1 \rightarrow E(Y_t)=\beta_1+(\beta_2+\beta_3)x_t\)
Considerato questo esempio la seguente applicazione lo generalizza fornendo un maggiore livello di dettaglio.
set.seed(1)
X <- runif(200,0, 15)
D <- sample(0:1, 200, replace = T)
Y <- 450 + 150 * X + 500*D +50 * (X * D) + rnorm(200, sd = 300)
m <- rbind(c(1, 2), c(3, 0))
graphics::layout(m)
#modello senza interazioni
plot(X, log(Y),
pch = 20,
col = "steelblue",
main = "Different Intercepts, Same Slope")
mod1_coef <- lm(log(Y) ~ X + D)$coefficients
abline(coef = c(mod1_coef[1], mod1_coef[2]),
col = "red",
lwd = 1.5)
abline(coef = c(mod1_coef[1] + mod1_coef[3], mod1_coef[2]),
col = "green",
lwd = 1.5)
#modello base più interazioni
plot(X, log(Y),
pch = 20,
col = "steelblue",
main = "Different Intercepts, Different Slopes")
mod2_coef <- lm(log(Y) ~ X + D + X:D)$coefficients
abline(coef = c(mod2_coef[1], mod2_coef[2]),
col = "red",
lwd = 1.5)
abline(coef = c(mod2_coef[1] + mod2_coef[3], mod2_coef[2] + mod2_coef[4]),
col = "green",
lwd = 1.5)
# dummy solo nell'interazione
plot(X, log(Y),
pch = 20,
col = "steelblue",
main = "Same Intercept, Different Slopes")
mod3_coef <- lm(log(Y) ~ X + X:D)$coefficients
abline(coef = c(mod3_coef[1], mod3_coef[2]),
col = "red",
lwd = 1.5)
abline(coef = c(mod3_coef[1], mod3_coef[2] + mod3_coef[3]),
col = "green",
lwd = 1.5)
Le interazioni si possono presentare in 4 forme:
Se \(\beta_3\) non è significativamente diverso da 0, allora l’interazione tra \(x_1\) e \(x_2\) non è utile per spiegare la variazione della variabile dipendente. In questo caso, la relazione tra \(x_1\) e \(x_2\) è detta additiva;
Se \(\beta_3\) è minore di 0 e indipendentemente le variabili \(x_1\) e \(x_2\) influenzano la variabile dipendente, allora la relazione è detta antagonistica;
Se \(\beta_3\) è maggiore di 0 e indipendentemente le variabili \(x_1\) e \(x_2\) influenzano la variabile dipendente, allora la relazione è detta sinergica;
Se \(\beta_3\) è significativamente diverso da 0 ma \(x_1\) e/o \(x_2\) non influenzano la variabile dipendente, allora la relazione è detta atipica.
La feature engineering si pone come obiettivo la creazione di specifiche variabili esplicative che permettono di migliorare l’efficacia di un modello utilizzando le informazioni predittive più rilevanti, cercando di individuare le interazioni principali tra i regressori.
La conoscenza del contesto di studio è un punto fondamentale per guidare il processo di selezione dei termini di interazione, se questa non è disponibile risulta fondamentale definire principi che guidino il processo di selezione delle interazioni, costruendo un framework che permetta di definire le relazioni di causalità tra i predittori e la variabile dipendente. Wu e Hamada forniscono un framework per identificare le interazioni significative che si basa su 3 dimensioni principali:
La gerarchia di interazione: tanto più alto è il grado dell’interazione, quanto più bassa è la capacità di quest’ultima di spiegare le variazioni della variabile dipendente. E’ consigliabile considerare inizialmente interazioni formate da coppie di variabili per poi aumentare il grado di interazione. Questo permette anche una maggiore interpretabilità del modello.
Effetto sparsità: solo una parte dei possibili termini di interazione è utile a spiegare una quantità significativa della variazione della variabile dipendente.
Principio di ereditarietà: un’interazione può essere valutata solo se i termini di ordine inferiore che la compongono sono efficaci nello spiegare la variazione della variabile dipendente. Ci sono due tipologie: ereditarietà forte (tutti i termini dell’interazione devono essere significativi) ed ereditarietà debole (almeno un termine dell’interazione deve essere significativo).
Prima di ricercare le possibili interazioni tra regressori è necessario fare alcune considerazioni in modo da adottare la strategia migliore. Per prima cosa bisognerebbe capire se è possibile o meno valutare tutti i termini di interazione e scegliere il momento più opportuno per crearli; questo perché l’ordine delle operazioni di preprocessing influisce sulla capacità di trovare interazioni importanti.
library(caret)
library(glmnet)
library(tidymodels)
library(AmesHousing)
library(MASS)
library(gridExtra)
library(grid)
library(moments)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ggcorrplot)
ames <- make_ames()
set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)
Questo dataset contiene informazioni relative alla vendita di proprietà immobiliari in Ames(IA) dal 2006 al 2010.
Ad un primo approccio suscitano alcune domande come:
options(scipen=999)
ggplot(ames_train, aes(x=Sale_Price))+
geom_histogram(color="blue", fill="white", bins = 75)+
ggtitle("Qual è il costo di una casa?")
La proprietà più costosa è stata venduta ad un prezzo pari a 755000 Є mentre la casa meno costosa è stata venduta a 12789Є.
Il prezzo medio di vendità risulta pari a 180827.6Є e la mediana dei prezzi di vendità coincide con 162000Є.
La distribuzione presenta una leggera asimmetria positiva ed inoltre sono presenti numerosi outliers.
ggplot(ames_train, aes(x=Year_Built))+
geom_histogram(color="blue", fill="white", bins = 14)+
ggtitle("Quando sono state costruite le case?")
La maggior parte delle case sono state costruite dopo gli anni 50.
La casa più vecchia è stata costruita nel 1872.
La casa più nuova è stata costruita nel 2010 .
tbl_1 = ames_train_tbl %>%
dplyr::select(Year_Sold, Mo_Sold) %>%
group_by(Year_Sold, Mo_Sold ) %>%
count(Year_Sold, Mo_Sold) %>%
dplyr::mutate(V=paste( as.character(Year_Sold), " - ", as.character(Mo_Sold) ))
ggplot(tbl_1, aes( x = V , y=n ) ) +
geom_bar(stat="identity", fill="steelblue")+
theme(text = element_text(size=7),
axis.text.x = element_text(angle=90, hjust=1))+
xlab("") +
ggtitle("Quando sono state vendute le case?")
Emerge la presenza di un pattern legato alla stagionalità nelle vendite delle case, in particolare sono presenti dei picchi nei mesi di giugno e luglio.
tbl_2 = ames_train_tbl %>%
dplyr::select(Neighborhood) %>%
arrange(Neighborhood) %>%
group_by(Neighborhood) %>%
count()
ggplot(tbl_2, aes( x = Neighborhood , y=n ) ) +
geom_bar(stat="identity", fill="steelblue")+
theme(text = element_text(size=7),
axis.text.x = element_text(angle=90, hjust=1))+
xlab("") +
ggtitle("In che quartiere sono situate le case?")
La maggior parte delle case sono situate nei quartieri di North Aes, College_Creek e Old Town, mentre i quartieri di Greens, Green_hills e Landmark non presentano un gran numero di case.
ggplot(ames_train, aes(x=Gr_Liv_Area))+
geom_histogram(color="blue", fill="white", bins = 21)+
ggtitle("Quanto sono grandi le case? (in sq feet)")
In media una casa è grande 1500.5 sq ft, la mediana della grandezza delle case risulta pari a 1452 sq ft.
La casa più grande è di 5642 sq ft, mentre la più piccola 334 sq ft.
Esaminiamo la distribuzione della variabile dipendente, SalePrice.
ggplot(ames_train_tbl, aes(x=Sale_Price))+
geom_histogram(aes(y=..density..), bins = 75, colour="black", fill="white")+
geom_density(alpha=.2, fill="#FF6666")
La distribuzione di SalePrice è caratterizzata dalla presenza di una leggera asimetria positiva. Infatti, si ha:
skewness(ames_train$Sale_Price)
## [1] 1.661015
kurtosis(ames_train$Sale_Price)
## [1] 7.657046
In fase di analisi verranno prese in considerazione solo le seguenti variabili:
ggcorrplot(cor(num_features_tbl), hc.order = TRUE, type = "lower",
lab = TRUE)
Dalla matrice di correlazione non risultano variabili particolarmente correlate.
Per dataset che hanno un numero piccolo o moderato di predittori, è possibile valutare tutte le coppie di interazioni.
Tuttavia maggiore è il numero di termini di interazione che vengono valutati, maggiore è la probabilità di trovare delle interazioni dette false positive, ovvero che risultano statisticamente significative ma non a causa di una vera relazione e possono dunque peggiorare le performance predittive del modello.
Gli approcci tradizionali per individuare i termini di interazioni più rilevanti si basano sui modelli nidificati.
Si consideri un modello lineare con 2 soli regressori, si definisce main effects model il seguente modello:
\(M_1: y=\beta_0+\beta_1x_1+\beta_2x_2+\epsilon\)
E’ possibile costruire un secondo modello che aggiunge una potenziale interazione tra i due regressori.
\(M_2: y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2+\epsilon\)
Questi due modelli sono detti “nidificati” perché il primo modello è incluso nel secondo (\(M_1 \subseteq M_2\)).
In presenza di questa struttura, è possibile effettuare un confronto statistico sulla quantità di informazione addizionale che è stata acquisita dal termine di interazione.
In questo contesto, con riferimento a una regressione lineare, è necessario confrontare l’errore residuo di \(M_1\) e \(M_2\) per valutare se il miglioramento dell’errore,ponderato per i gradi di libertà, è sufficiente per essere considerato reale.
Per la regressione lineare, la funzione obiettivo utilizzata per confrontare i modelli è la verosimiglianza statistica (in questo caso, l’errore residuo).
Per altri modelli, ad esempio il modello di regressione logistica, la funzione obiettivo per confrontare i modelli nidificati coincide con la verosimiglianza binomiale.
Il risultato di questo test statistico lo si legge dal p-value che può essere interpretato come il tasso di risultati falsi positivi e riflette la probabilità che l’informazione aggiuntiva catturata dal termine di interazione sia dovuta alla casualità.
Più è basso il p-value minore è la probabilità che l’informazione aggiuntiva catturata dal termine di interazione sia dovuta alla casualità.
Una evoluzione del modello nidificato sfrutta la cross-validation in modo da creare più coppie di modelli annidati a partire da differenti versioni del training set, a differenza dell’approccio tradizionale in cui gli stessi dati usati per la creazione del modello vengono riutilizzati per la sua valutazione che può essere svolta mediante qualsiasi misura di performance.
Considerato che la valutazione dipende molto dal contesto, questi due approcci non garantiscono conformità nei risultati.
Come detto in precedenza più sono i confronti che vengono effettuati, più alta è la possibilità di trovare interazioni false-positive.
Esistono diversi metodi per affrontare questo problema, ad un estremo, si sceglie di non effettuare controlli per i falsi positivi, oppure si applica la correzione di Bonferroni che utilizza una “severa” penalità esponenziale per minimizzare eventuali risultati falsi positivi.
Tuttavia la tecnica più equilibrata per individuare e gestire i falsi positivi si basa sul false discovery rate (FDR).
Nella presente applicazione vengono considerati solamente 18 regressori presenti nel dataset Ames, rispettivamente: Bldg_Type, Neighborhood, Year_Built, Gr_Liv_Area, Full_Bath, Year_Sold, Lot_Area, Central_Air, Longitude, Latitude, MS_SubClass, Alley, Lot_Frontage, Pool_Area, Garage_Finish, Foundation, Land_Contour e Roof_Style.
Il comando recipe() permette di applicare opportune trasformazioni ai dati iniziali.
ames_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_bs(Longitude, Latitude, options = list(df = 5)) %>%
step_zv(all_predictors())
In particolare, le trasformazioni effettuate sono state le seguenti:
## Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 18
##
## Operations:
##
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood
## B-Splines on Longitude, Latitude
## Zero variance filter on all_predictors()
Si procede creando una dataframe contenente tutte le possibile coppie di interazioni \(\frac{p(p-1)}{2}\) tra i regressori.
ames_preds <- ames_rec$var_info$variable[ames_rec$var_info$role == "predictor"]
interactions <- t(combn(ames_preds, 2))
colnames(interactions) <- c("var1", "var2")
Si decide di addestrare un modello lineare sui dati del training set utilizzando una repeated k-fold cross-validation, iterando la procedura 5 volte con K=10.
int_ctrl <- trainControl(method = "repeatedcv", repeats = 5)
set.seed(2691)
main_eff <- train(ames_rec,
data = ames_train,
method = "lm",
metric = "RMSE",
trControl = int_ctrl)
Si implementa la funzione compare_models_1way al fine di confrontare il modello senza interazioni (main effect model) e il main effect model più un’ interazione.
compare_models_1way <- function(a, b, metric = a$metric[1], ...) {
mods <- list(a, b)
rs <- resamples(mods)
diffs <- diff(rs, metric = metric[1], ...)
diffs$statistics[[1]][[1]]
}
In questo modo è possibile valutare l’effetto delle interazioni in termini di: Reduction, Pvalue, RMSE, anova_p.
for (i in 1:nrow(interactions)) {
tmp_vars <- c("Class", interactions$var1[i], interactions$var2[i])
int <-
as.formula(paste0(
"~ starts_with('",
interactions$var1[i],
"'):starts_with('",
interactions$var2[i],
"')"
))
int_rec <-
ames_rec %>%
step_interact(int) %>%
step_zv(all_predictors())
set.seed(2691)
main_int <- train(int_rec,
data = ames_train,
method = "lm",
metric = "RMSE",
trControl = int_ctrl)
tmp_diff <- compare_models_1way(main_int, main_eff, alternative = "less")
interactions$RMSE[i] <- getTrainPerf(main_eff)[1, "TrainRMSE"]
interactions$Reduction[i] <- -tmp_diff$estimate
interactions$Pvalue[i] <- tmp_diff$p.value
a1 <-
anova(main_eff$finalModel, main_int$finalModel) %>%
tidy() %>%
slice(2) %>%
pull(p.value)
interactions$anova_p[i] <- a1
}
head(interactions)
## # A tibble: 6 × 6
## var1 var2 Reduction Pvalue RMSE anova_p
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Bldg_Type Neighborhood 0.000184 0.0479 0.0748 2.08e- 4
## 2 Bldg_Type Year_Built -0.0000263 0.680 0.0748 7.94e- 2
## 3 Bldg_Type Gr_Liv_Area 0.000736 0.000000471 0.0748 7.14e-10
## 4 Bldg_Type Full_Bath 0.000215 0.00107 0.0748 1.53e- 3
## 5 Bldg_Type Year_Sold -0.0000560 0.915 0.0748 2.33e- 1
## 6 Bldg_Type Lot_Area -0.0000249 0.683 0.0748 6.88e- 2
plot <-
ggplot(all_p, aes(x = Value)) +
geom_histogram(breaks = (0:20)/20) +
geom_rug(col = "red") +
facet_grid(Estimate ~ Method)
Questa rappresentazione confronta la distrubuzione dei p-value ottenuta con l’approccio tradizionale e la cross-validation, valutando le diverse tipologie di aggiustamento (No Adjustment, FDR, Bonferroni).
Emerge che utilizzando l’approccio tradizionale quasi tutti gli effetti delle interazioni vengono valutate significative, ciò potrebbe essere dovuto a una potenziale presenza di overfitting. L’approccio di cross validation ha un tasso minore di significatività ed è meno influenzato da entrambi i tipi di aggiustamenti.
Valutare i termini di interazione singolarmente impedisce ai modelli semplici di valutare l’importanza dei termini di interazione nel modello completo ( contenente tutti i main effects e i relativi termini di interazione), per questo se la completa enumerazione è possibile si potrebbe provare un approccio che ricerca le interazioni partendo da un modello completo.
Utilizzando quest’approccio potrebbe però capitare che nei dati si presentino più predittori rispetto alle osservazioni del campione, ciò potrebbe non essere un problema per modelli come basati su alberi decisionali, neural networks, SVM, e K-nearest neighbors, che sono applicabili anche quando i dati contengono più predittori rispetto alla numerosità delle osservazioni del campione. Tuttavia altre tecniche più interpretabili e che spesso riportano buone performance, ad esempio le regressioni lineari e logistiche, non possono essere utilizzate direttamente sotto queste condizioni. Questa limitazione deriva dal fatto che qualora si hanno più predittori che osservazioni nel campione, oppure un predittore può essere scritto come combinazione di altri, non è possibile effettuare l’inversione della matrice del disegno.
Essendo l’interpretabilità e la bontà del modello obiettivi della feature engineering, si sono cercati metodi per utilizzare anche in questo caso la regressione lineare e logistica. In questo contesto, si introduce una famiglia dei tecniche di modellazione dette modelli penalizzati.
I metodi più utilizzati tra i modelli penalizzati sono la ridge regression e la lasso regression. Questi due tipi differenti di penalità sono solitamente associati a differenti task:
Talvolta potrebbe essere necessario combinare le due penalità, a tal fine è stato introdotto il modello glmnet.
Nel modello glmnet compaiono entrambi i termini di penalizzazione ridge e lasso, inoltre, sono presenti due parametri di tuning:
Ad esempio scegliendo \(\alpha=1\) avremo un modello fully lasso penality, scegliendo \(\alpha=0.5\) avremo un modello che è metà lasso penality e metà ridge penality.
La regressione tramite modello glmnet si può vedere, come nel caso della regressione lineare, come un problema di minimizzazione, in particolare del valore del somma dei quadrati dei residui, calcolata con la seguente equazione:
\[ SSE=\sum_{i=1}^n (y_i - \hat{y}_i)^2 +(1-\alpha)\lambda_r \sum_{j=1}^P \beta^2_j + \alpha \lambda_l\sum_{j=1}^P |\beta_j| \]
L’esempio è stato svolto utilizzando i dati ames e consiste nel tuning di due modelli glmnet. Un primo modello contenente solo i main effect e un secondo modello contenente main effect e relative interazioni. Il modello glmnet dopo aver provato tutte le combinazioni di valori dei parametri di tuning \(\alpha\) e \(\lambda\) sceglie il “best tuning”, cioè il modello che minmizza il valore dell’RMSE. La scelta dei parametri di tuning ottimali permette di selezionare le variabili esplicative più significative e di stimare i loro coefficienti.
Il primo passo è quello di recuperare i dati ames e di effettuare lo split in training set e test set e successivamente con il comando vfold_cv() in 10 folds per poi svolgere una cross validation.
Si passa successivamente alla creazione del primo modello, come detto in precedenza, creato utilizzando solo i main effect tramite una struttura chiamata recipe() che permette inoltre di specificare una pipeline per il preprocessing dei dati. In questo passo si specificano quindi la variabile target e quelle esplicative del modello e le operazioni da effettuare in fase di preprocessing.
main_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_zv(all_predictors()) %>%
step_bs(Longitude, Latitude, options = list(df = 5)) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors())
Dopo aver creato un dataframe con tutte le possibili interazioni di secondo ordine tra i main effects del primo modello, si crea una nuova recipe() per il secondo modello.
int_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_interact(interactions) %>% # Aggiunta delle interazioni
step_zv(all_predictors()) %>%
step_bs(Longitude, Latitude, options = list(df = 5)) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors())
Il modello glmnet, utilizza due diversi parametri di tuning \(\alpha\) e \(\lambda\), si crea un griglia contenente tutte le possibili combinazioni di questi parametri che il modello utilizzerà per effettuare il tuning ricercando la coppia che minimizza il valore dell’RMSE. Il parametro \(\alpha\) può assumere valori da \(0.2\) a \(1\), con passo di \(0.2\). Mentre \(\lambda\) valori da \(10^{-4}\) a \(10^{-1}\) con il valore dell’esponente che varia di \(0.1\).
glmn_grid <- expand.grid(alpha = seq(.2, 1, by = .2), lambda = 10^seq(-4, -1, by = 0.1))
main_glmn <-
train(main_rec, # La recipe
data = ames_train, # Dati da utilizzare per il tuning dei parametri
method = "glmnet",
tuneGrid = glmn_grid, # Dataframe con le combinazioni dei parametri di tuning
trControl = ctrl # Divisione in fold per la cross validation
)
## $all
## [1] 58
##
## $main
## [1] 42
##
## $perf
## TrainRMSE TrainRsquared TrainMAE method
## 1 0.0754791 0.8184643 0.05187811 glmnet
Dopo aver effettuato il tuning dei parametri per il primo modello, il modello risultante che minimizza il valore dell’RMSE è composto da 42 main effects che sono stati selezionati tramite la penalizzazione lasso.
int_glmn <-
train(int_rec,
data = ames_train,
method = "glmnet",
tuneGrid = glmn_grid,
trControl = ctrl
)
## $all
## [1] 981
##
## $main
## [1] 9
##
## $int
## [1] 132
##
## $perf
## TrainRMSE TrainRsquared TrainMAE method
## 1 0.07589794 0.817775 0.05157754 glmnet
Il secondo modello ottimale selezionato tramite i parametri di tuning è composto invece da 9 main effects e 132 termini di interazione.
Vengono mostrati i parametri che permettono di minimizzare il valore dell’RSME.
int_glmn$bestTune
## alpha lambda
## 50 0.4 0.006309573
Il plot che mostra il rapporto tra il valore dell’RMSE (asse y) e quello di \(\alpha\) (colori) e \(\lambda\) (asse x).
Dopo aver effettuato il tuning dei precedenti due modelli e trovato per ogni modello i regressori significativi, l’esempio procede con un modello glmnet two way. Per la sua creazione vengono presi i main effects risultati significativi nel primo modello, si creano tutte le possibili interazioni binarie e infine si aggiungono ad un modello contenente i main effect. Si riesegue poi il tuning su questo terzo modello.
two_stage_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_interact(interaction_subset) %>%
step_zv(all_predictors()) %>%
step_bs(Longitude, Latitude, options = list(df = 5)) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors())
## $all
## [1] 741
##
## $main
## [1] 9
##
## $int
## [1] 108
##
## $perf
## TrainRMSE TrainRsquared TrainMAE method
## 1 0.07614744 0.8166224 0.05192411 glmnet
Il numero di interazioni da studiare cresce esponenzialmente al crescere del numero di variabili esplicative rendendo l’esplorazione dei termini di interazione molto onerosa dal punto di vista computazionale.
Dal lato statistico non sarebbe una scelta saggia studiare tutte le possibili interazioni. Le interazioni importanti tra le variabili esplicative, infatti, sono generalmente poche e all’aumentare del numero di interazioni irrilevanti esaminate aumenta la penalità imposta al fine di valutarle. Di conseguenza, interazioni realmente significative potrebbero essere perse. In questi casi sono necessari altri approcci per rendere la ricerca più efficiente, ma comunque efficace, nel trovare le interazioni importanti.
I metodi presentati finora sono efficaci per individuare interazioni in presenza di un numero moderato di regressori. Metodi basati sugli alberi decisionali risultano particolarmente efficaci per scoprire potenziali interazioni in presenza di molti regressori.
Queste tecniche si basano sul partizionamento ricorsivo che consiste nell’individuare il punto di split ottimo di un predittore che separa le osservazioni in gruppi omogenei rispetto alla variabile risposta.
Una volta effettuato lo split, la stessa procedura si ripete per ciascun sottoinsieme di osservazioni in corrispondenza di ogni nodo successivo.
Questo processo termina quando si raggiunge il numero minimo di osservazioni in ciascun nodo oppure viene soddisfatto un criterio di ottimizzazione.
Quando i dati vengono splittati ricorsivamente, ogni nodo successivo può essere interpretato come una interazione locale tra il nodo precedente e il nodo corrente.
Inoltre, più è frequente il verificarsi di nodi successivi della stessa coppia di predittori, maggiore è la probabilità che l’interazione tra questi sia di tipo globale.
Formalmente, al fine di comprendere l’efficacia della tecnica presentata per individuare una relazione sinergica tra due regressori, si consideri un esempio in cui l’albero di partizionamento ricorsivo ottimo (Figura 1) risulta:
Figura 1: Albero di partizione ricorsiva utile ad individuare una relazione sinergetica tra due predittori
L’interazione viene catturata attraverso diversi splits alternati tra i due predittori ai livelli successivi dell’albero.
Considerando il loro tipo di equazioni, i modelli basati sugli alberi decisionali possono essere interpretati come una interazione pura.
Ad esempio, l’equazione del nodo 4 può essere scritta come:
\[node_4= I(x_1 <0.655)*I(x_2<0.6)*I(x_0.316)*0.939 \]
dove 0.939 coincide con la media delle osservazioni presenti nel nodo 4 mentre \(I\) rappresenta la funzione indicatrice.
Una relazione sinergica produce delle curve di livello non lineari per la variabile risposta (Figura 2), il modello di partizionamento ricorsivo rompe lo spazio in regioni rettangolari, pertanto risultano necessarie più regioni al fine di spiegare correttamente la risposta, questo però limita fortemente la capacità di modellare interazioni globali.
Figura 2 : Risposta prevista ottenibile a partire da un modello di partizione ricorsiva tra 2 regressori che sono in relazione sinergica
Questo grafico rappresenta le previsioni di ciascun nodo terminale rappresentate da regioni rettangolari in cui il colore indica il valore della risposta prevista per la corrispondente regione a cui sono sovrapposte le curve di livello generate dal modello lineare utilizzato per stimare il termine di interazione.
Nonostante il modello stimato, decomponendo lo spazio della variabili in regioni rettangolari, approssimi adeguatamente la relazione sinergica, sono richieste troppe regioni per spiegare questa relazione.
Un modello basato su un singolo albero decisionale può essere caratterizzato da un’alta variabilità, ciò rende le previsioni poco affidabili.
Per superare questo limite, è possibile ricorrere a metodi di ensemble learning come bagging e boosting che riescono ad identificare le principali relazioni tra i predittori, garantendo una migliore capacità previsiva della variabile risposta.
Le tecniche di ensemble learning sono particolarmente efficaci nell’identificare le relazioni tra predittori e risposta perché aggregano molti alberi decisionali utilizzando delle versioni leggermente diverse dei dati originali.
Considerata una relazione sinergica tra x1 e x2, la Figura 3 mostra che i modelli boosted tree e bagged tree sono in grado di approssimare meglio la relazione tra i predittori e la variabile risposta rispetto a un singolo albero decisionale.
Figura 3 : Modello boosted tree vs. Modello bagged tree nel caso di una relazione sinergica tra due predittori
Infatti, la radice quadrata dell’errore quadratico medio (RMSE) per questi 3 modelli risulta:
\(RMSE_{boosting}=0.27\)
\(RMSE_{bagging}=0.6\)
\(RMSE_{tree}=0.9\)
Nonostante la loro complessità e le loro difficoltà nell’identificare le interazioni globali, quando si verificano interazioni importanti all’interno di uno o più sottoinsiemi dei campioni i metodi basati sugli alberi decisionali risultano particolarmente adeguati per modellare questo tipo di interazioni (interazioni locali).
In conclusione, gli alberi decisionali e i metodi di ensemble learning permettono di ottenere delle informazioni rilevanti riguardo le principali interazioni tra regressori, è suggeribile utilizzare tali informazioni in modo da costruire un modello lineare più semplice.
Un altro approccio basato sui metodi di ensemble learning introdotto da Friedman & Popescu (2008) utilizza la dipendenza parziale al fine di identificare le principali interazioni.
In altre parole, questa tecnica confronta l’effetto congiunto di due (o più) predittori con l’effetto individuale di ciascun predittore di un modello.
Se un predittore non interagisce con nessun altro allora la differenza tra l’effetto congiunto e quello individuale sarà vicina a 0.
Qualora, invece, fosse presente un interazione tra un predittore ed un altro allora questa differenza sarebbe maggiore di 0.
Questo confronto avviene per mezzo della statistica H.
Questa statistica presenta distribuzione bimodale qualora si è in presenza di valori sparsi, e solitamente, ciò accade quando un predittore non è coinvolto in nessuna interazione.
Quando il dataset a disposizione è relativamente piccolo, si calcola la statistica H come la mediana delle statistiche H relative a più campioni di tipo bootstrap.
Successivamente viene fissato un cut off relativamente alla statistica H in modo da selezionare le variabili.
Le variabili a cui corrisponde un valore della statistica H maggiore del cut off vengono segnalate come coinvolte in una potenziale interazione.
In questo modo, è possibile inserire in un modello lineare tutte le coppie di interazioni segnalate.
Per comprendere la metodologia presentata si consideri la seguente applicazione:
h_stat <- function(data, ...) {
require(foreach)
pre_obj <- pre(data = data, ...)
pred_names <- pre_obj$x_names
pred_df <-
tibble(
Predictor = pred_names,
H = 0
)
int_obj <-
try(
interact(
pre_obj,
parallel = TRUE,
plot = FALSE
),
silent = TRUE)
if (!inherits(int_obj, "try-error")) {
res <-
tibble(
Predictor = names(int_obj),
H = as.numeric(int_obj)
)
missing_pred <- setdiff(pred_names, res$Predictor)
if (length(missing_pred) > 0) {
res <-
pred_df %>%
dplyr::filter(Predictor %in% missing_pred) %>%
bind_rows(res)
}
} else {
print(int_obj)
res <- pred_df
}
res
}
boot_h <- function(split, ...) {
dat <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = analysis(split)) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, MS_SubClass, Roof_Style, Foundation,
threshold = 0.05) %>%
step_zv(all_predictors()) %>%
prep(training = analysis(split), retain = TRUE) %>%
juice()
h_stat(data = dat, ...)
}
Per motivi computazionali, in questa applicazione vengono considerati solamente 18 regressori presenti nel dataset Ames, rispettivamente: Bldg_Type, Neighborhood, Year_Built, Gr_Liv_Area, Full_Bath, Year_Sold, Lot_Area, Central_Air, Longitude, Latitude, MS_SubClass, Alley, Lot_Frontage, Pool_Area, Garage_Finish, Foundation, Land_Contour e Roof_Style.
ames <- make_ames()
set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)
ames_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, MS_SubClass, Roof_Style, Foundation,
threshold = 0.05) %>%
step_zv(all_predictors())
ames_pre <- ames_rec %>% prep(ames_train) %>% juice()
ames_rec
## Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 18
##
## Operations:
##
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood, MS_SubClass, Roof_Style, Foundation
## Zero variance filter on all_predictors()
Il comando recipe() permette di applicare opportune trasformazioni ai dati iniziali.
Dopo aver trasformato le variabili si calcola la statistica H sul training set.
set.seed(3108)
obs_h <-
h_stat(
data = ames_pre,
Sale_Price ~ .,
normalize = TRUE,
ntrees = 10,
winsfrac = 0,
verbose = TRUE,
par.final = TRUE
)
Considerato che il dataset a disposizione è relativamente piccolo, si procede calcolando la statistica H anche come la mediana delle statistiche H relative a più campioni di tipo bootstrap.
h_stats <-
h_vals %>%
pull(stats) %>%
bind_rows() %>%
group_by(Predictor) %>%
summarise(Bootstrap = median(H)) %>%
full_join(
obs_h %>% dplyr::rename(`Original Data` = H)
) %>%
arrange(Bootstrap)
ordered_pred <-
h_stats %>%
pull(Predictor)
h_stats <-
h_stats %>%
gather(Estimator, H, -Predictor) %>%
mutate(Predictor = factor(Predictor, levels = ordered_pred)) %>%
dplyr::filter(H > 0)
Successivamente, viene fissato un cut off pari a 0.001 in modo da segnalare le variabili che risultano coinvolte in una potenziale interazione.
p <-
h_stats %>%
mutate(
Predictor = clean_value(Predictor),
Predictor = str_replace(Predictor, ":", "")
) %>%
ggplot(aes(x = reorder(Predictor, H), y = H, shape = Estimator)) +
geom_hline(yintercept = 0.001, col = "red", alpha = 0.5, lty = 2) +
geom_point() +
coord_flip() +
xlab("") +
scale_shape_manual(values = c(16, 1)) +
theme(legend.position = "top")
p
Da questa rappresentazione grafica emerge che 6 predittori sono coinvolti in almeno una interazione.
Un secondo approccio utile ad identificare le interazioni combina le informazioni riguardo l’importanza dei predittori con i principi di gerarchia, sparsità ed ereditarietà. Se una interazione viene considerata importante, allora un metodo basato sugli alberi di decisione utilizzerà i predittori coinvolti più volte in più alberi in modo da scoprire la loro relazione con la variabile risposta.
I coefficienti di regressione associati a tali predittori risulteranno significativamente diversi da 0.
In questo modo è possibile creare le principali coppie di interazioni e valutare la loro rilevanza nelle previsione della variabile risposta.
Ogni metodo basato sugli alberi decisionali è caratterizzato da uno specifico algoritmo utile a calcolare l’importanza di un predittore.
In questo contesto, uno degli algoritmi più popolari si ottiene a partire da un random forest che utilizza campioni bootstrap per combinare le previsioni ottenute da diversi modelli, riuscendo così a migliorare la capacità previsiva del modello.
Un albero decisionale costruito su un campione bootstrap utilizza \(\frac{2}{3}\) delle osservazioni originali.
E’ possibile prevedere la risposta della i-esima osservazione utilizzando gli alberi per i quali questa risulta Out-of-Bag.
Questo produce circa \(\frac{B}{3}\) previsioni per la i-esima osservazione che vengono aggregate calcolando la media oppure la moda.
Per valutare l’importanza di uno specifico predittore, i suoi valori vengono mescolati, e successivamente viene ricalcolato l’errore Out-of-Bag.
Qualora il predittore sia rilevante, tale errore dovrebbe peggiorare in modo significativo.
La random forest conduce questa analisi per ogni predittore presente nel modello e crea uno score per ciascuna feature, all’aumentare dello score aumenta il grado di importanza del predittore.
Innanzitutto, si procede stimando un modello di tipo random forest e si ottiene uno score che riflette l’importanza di ciascun regressore.
rf_mod <- ranger(Sale_Price ~ ., data = ames_pre, num.trees = 2000,
importance = "impurity", num.threads = workers,
seed = 8516)
rf_imp <-
tibble(
Predictor = names(rf_mod$variable.importance),
Importance = unname(rf_mod$variable.importance)
)
p <-
ggplot(rf_imp, aes(x = reorder(Predictor, Importance), y = Importance)) +
geom_bar(stat = "identity") +
coord_flip() +
xlab("")
p
Questa figura mostra i 10 regressori che risultano più importanti, alcuni di loro sono stati segnalati come potenziali termini di interazione da altre metodologie presentate.
I predittori più importanti sono living area e year built, è ragionevole assumure che una possibile interazione tra questi regressori risulti molto significativa.
In generale, le interazioni tra i principali predittori potrebbero essere incluse in un modello più semplice in modo da valutare se migliora la capacità previsiva della variabile risposta.
Il problema della selezione delle interazioni più importanti, anche considerato un numero moderato di features, risulta impossibile da affrontare con una valutazione completa di tutti i possibili modelli.
Quando si utilizzano modelli lineari si sceglie spesso di utilizzare metodi di forward, backward e stepwise selection.
Miller ha proposto una tecnica alternativa per risolvere i problemi del modello lineare. Anziché rimuovere un predittore, questa tecnica utilizza un approccio sostitutivo: ogni predittore selezionato è scelto sistematicamente per sostituirne un altro all’interno del modello.
Al fine di comprendere la tecnica presenta viene riportato un esempio in cui partendo da un modello formato da dieci variabili esplicative si cerca il miglior modello con solo tre predittori. Un’iterazione dell’algoritmo si può schematizzare come segue:
Non è garantito che l’ottimo trovato sia un ottimo globale, pertanto bisognerebbe ripetere la procedura con diverse partenze casuali, cioè con combinazioni di variabili iniziali diverse e determinare la soluzione ottima raggiunta.
Hawkins ha esteso questo algoritmo, creando un approccio di modellazione chiamato Feasible Solution Algorithm (FSA), uno dei principali vantaggi è la dimensione dello spazio di ricerca nell’ordine \(q*m*p\) , molto più piccolo dell’intero spazio \(p^m\).
L’algoritmo precedentemente presentato non considera le interazioni, Lambert lo ha generalizzato ampliando lo spazio di ricerca anche ad esse. L’approccio inizia con un modello base che include le variabili ritenute significative per la previsione della variabile target e, come prima cosa, viene scelto l’ordine delle interazioni tra i main effects. Il processo di identificazione delle interazioni segue la logica dell’FSA orginale, viene riportato il seguente esempio:
Le soluzioni/gli ottimi trovati con le diverse partenze casuali vengono utilizzati per identificare le potenziali soluzioni.
Anche per l’algoritmo FSA l’esempio viene effettuato sfruttando i dati ames. Per poter confrontare i risultati ottenuti con quelli delle procedure precedenti si è scelto di trasformare anche in questo caso le variabili nominali in variabili dummy individuali e poi procedere alla creazione delle interazioni (la scelta delle variabili dummy, per la creazione delle interazioni, è stata vincolata in modo che non possano essere accoppiate variabili derivanti dallo stesso predittore). Nell’implementazione dell’algoritmo sono state utilizzate 50 partenze casuali, con un subset di massimo 30 variabili esplicative, questo porta ad avere in ogni partenza casuale un massimo di 60 swaps delle variabili che compongono l’interazione.
Il primo passo consiste nel recuperare i dati ames ed effettuare lo split in training set e test set e successivamente con il comando vfold_cv() in 10 folds per poi svolgere una cross validation.
Tramite la struttura recipe() si crea la matrice delle variabili esplicative relativa al modello contenente solo main effect. Grazie allo stesso comando si crea la pipeline di preprocessing dei dati, in particolare le variabili nominali sono trasformate in variabili dummies individuali.
ames_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_bs(Longitude, Latitude, options = list(df = 5)) %>%
step_zv(all_predictors())
ames_rec
## Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 18
##
## Operations:
##
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood
## Dummy variables from all_nominal()
## B-Splines on Longitude, Latitude
## Zero variance filter on all_predictors()
Viene mostrato il contenuto della recipe(), in questo caso il modello è formato da una variabile target e 18 predittori.
Con il seguente comando viene creato un set di metriche che verranno poi utilizzate nell’esecuzione dell’algoritmo.
multi_metric <- metric_set(rmse, mae, rsq)
Un’altra informaznione necessaria per l’algoritmo è la specificazione del modello di regressione da utilizzare. Nel nostro caso è stata usata una regressione lineare.
lr_spec <- linear_reg() %>%
set_engine("lm")
Si esegue l’algoritmo fsa appoggiandosi alla funzione fsa_two_way contenuta nel file fsa_function. La funzione esegue un fsa two way, in quanto si cercano solo le interazioni di secondo grado tra i predittori. Per l’esecuzione della suddetta funzione vengono forniti i seguenti parametri:
la funzione ha come parametri di default:
set.seed(236)
ames_search <- fsa_two_way(ames_folds, ames_rec, lr_spec, multi_metric)
## # A tibble: 31 x 4
## var_1 var_2 pval RMSE
## <chr> <chr> <dbl> <dbl>
## 1 Foundation: CBlock Living Area 0.00539 0.0513
## 2 Foundation: PConc Living Area 0.00828 0.0515
## 3 Building Type Duplex Living Area 0.00178 0.0516
## 4 Living Area MS SubClass: Duplex All Styles and Ages 0.00102 0.0516
## 5 Roof Style: Hip Year Built 0.0495 0.0518
## 6 Roof Style: Gable Year Built 0.000320 0.0518
## 7 Longitude (spline) Neighborhood: Northridge Heights 0.00789 0.0518
## 8 Living Area Neighborhood: Northridge Heights 0.0964 0.0519
## 9 Full Bath Land Contour: HLS 0.0767 0.0519
## 10 Foundation: CBlock Lot Area 0.0170 0.0519
## # ... with 21 more rows
L’output dell’algoritmo mostra le interazioni ordinate per il valore dell’RSME ottenuto dal modello con i main effect più l’interazione in questione. Inoltre è presente anche il valore del p-value ottenuto nel t-test per verificare se il miglioramento apportato dall’aggiunta dell’interazione al modello base si può ritenere significativo.
Per esempio sembra esserci una potenziale interazione tra il building type e la living area (p-value basso 0,00178).
Esistono altri modelli utili a scoprire le principali interazioni tra regressori basati sulla Multivariate adaptive regression splines (MARS) e sul modello Cubist.
In generale la Multivariate adaptive regression splines (MARS) è un modello di regressione non lineare per una variabile risposta di tipo continuo, essa ha come obiettivo la creazione di una hinge function utile a descrivere la relazione tra il predittore e la variabile dipendente. Tale funzione consente al modello di cercare relazioni non lineari. Nel contesto di ricerca dei principali termini di interazione, la MARS può essere utilizzata per cercare dei prodotti tra regressori utili a creare delle interazioni non lineari che isolano porzioni dello spazio delle variabili.
Cubist è un modello di regressione rule-based che costruisce un albero iniziale e lo scompone in un insieme di regole che vengono potate e qualora necessario eliminate. Per ciascuna regola viene definito un apposito modello lineare, creando un insieme di interazioni disgiunte. Le previsioni possono appartenere a più regole e, in questo caso, viene calcolata la media delle previsioni del modello lineare associato. Questa struttura è molto flessibile e ha una forte probabilità di individuare le interazioni locali all’interno dell’intero spazio dei predittori.
L’errore di previsone \(ErrF\) è definito come:
\[ErrF=E(MSE_{Te})=E \left[\frac{1}{n}\sum_{i=1}^{n} (Y^*_i-\hat f(x_i))^2\right] = \frac{1}{n}\sum_{i=1}^{n} E\left [(Y^*_i-\hat f(x_i))^2\right]\] E’ possibile distinguere tre tipologie di fonti che causano un errore di previsione:
Errore riducibile, scomponibile in:
Distorsione: esprime la quota di errore legata alla differenza media tra \(f\) e \(\hat f\)
Varianza: variabilità dello stimatore \(\hat f\)
Non esiste un modello che minimizza contemporaneamente la distorsione e la varianza, è necessario cercare un “compromesso”.
\[E[(Y^*_i-\hat Y^*_i)^2]=E[(f(x_i) + \epsilon^*_i - \hat f(x_i))^2]= \underbrace{E[(f(x_i) - \hat f(x_i))^2]}_{riducibile}+\underbrace{Var(\epsilon^*_i)}_{non\ riducibile} \] dove, \(Var(\epsilon^*_i)=\sigma^2\).
L’errore riducibile può essere scomposto nel \([Bias(\hat f(x_i))]^2\) e nella varianza \(Var( \hat f(x_i))\) dello stimatore \(\hat f\), rispettivamente:
\[E[ ( f(x_i) - \hat f(x_i) )^2 ] = \underbrace{[ E ( \hat f(x_i) - f(x_i) ) ]^2}_{[Bias( \hat f(x_i) ) ]^2 } + \underbrace{Var[\hat f(x_i)]}_{Variance[\hat f(x_i)]}\]
A partire da questa scomposizione l’errore di previsione \(ErrF\) si può scrivere come:
\[ErrF=\sigma^2+\underbrace{\frac{1}{n}\sum_{i=1}^{n} (E\hat f(x_i)-f(x_i))^2}_{Bias^2}+\underbrace{\frac{1}{n} \sum_{i=i}^{n} Var(\hat f(x_i) ) }_{Varianza}\]
Bias e varianza non possono essere minimizzate simultaneamente, qualora si fosse interessati a diminuire il bias si dovrà utilizzare un modello molto complesso ma in questo modo si otterrà un aumento della varianza. Esiste un trade-off tra bias e varianza:
Considerare il dataset ames come Population_data.
Dividere Population_data in Test_data e Training_data.
Costruire Population_model: addestrato su Population_data e applicato su Test_data.
Costruire il Mean_model, a partire da 20 campioni casuali estratti da Training_data è stato costruito un modello per ognuno di essi. Ciò permette di ottenere la previsione media dei modelli sul Test_data.
Calcolo Bias
Calcolo Varianza
library(caret)
library(glmnet)
library(tidymodels)
library(AmesHousing)
library(MASS)
library(ipred)
library(rpart)
library(randomForest)
library(reshape2)
\[Variance (x) = E_{sample}[\hat f_{sample}(x) -\bar f(x)]^2\] Dove,
\(\bar f(x)\) modello medio, media delle previsioni ottenute a partire da 20 modelli costruiti su 20 campioni differenti.
\(\hat f_{sample}(x)\) valori previsti sul test set, stimati a partire da un modello costruito su un campione.
calculate_variance_of_model <- function(samplePredictions, y_test){
predictions_mean_model <- colMeans(samplePredictions)
colNames <- colnames(samplePredictions)
variance = numeric()
mse = numeric()
i = 1
for (colName in colNames) {
variance[i] = mean(as.numeric(samplePredictions[,colName] - predictions_mean_model[i])^2)
mse[i] = mean((samplePredictions[, colName] - as.numeric(y_test[i,]))^2)
i=i+1
}
return(list(mean(variance), mean(mse)))
}
\[Bias (x) = (f(x) -\bar f(x))\]
Dove,
\(\bar f(x)\) modello medio, media delle previsioni ottenute a partire da 20 modelli costruiti su 20 campioni differenti.
\(f(x)\) modello generatore dei dati.
calculate_bias_of_model <- function(samplePredictions, y_hat_pop){
predictions_mean_model <- colMeans(samplePredictions)
return((mean(abs(predictions_mean_model-y_hat_pop)))^2)
}
ames <- make_ames()
set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)
load("C:/Users/Stefano/Downloads/ames_h_stats (1).RData")
int_vars <-
h_stats %>%
dplyr::filter(Estimator == "Bootstrap" & H > 0.001) %>%
pull(Predictor)
interactions <- t(combn(as.character(int_vars), 2))
colnames(interactions) <- c("var1", "var2")
interactions <-
interactions %>%
as_tibble() %>%
mutate(
term =
paste0(
"starts_with('",
var1,
"'):starts_with('",
var2,
"')"
)
) %>%
pull(term) %>%
paste(collapse = "+")
interactions <- paste("~", interactions)
interactions <- as.formula(interactions)
int_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_interact(interactions) %>%
step_zv(all_predictors()) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors()) %>%
step_bs(Longitude, Latitude, options = list(df = 5))
prep.int_rec <- prep(int_rec)
train_data <- bake(prep.int_rec, new_data = ames_train)
test_data <- bake(prep.int_rec, new_data = ames_test)
ames_data <- bake(prep.int_rec, new_data = ames)
#x_train <- train_data[,-8]
#y_train <- train_data[,8]
x_test <- test_data[,-8]
y_test <- test_data[,8]
#x_pop <- ames_data[,-8]
#y_pop <- ames_data[,8]
Considerato che la funzione generatrice dei dati è ignota, si assume che tale coincida con il modello stimato sulla totalità dei dati a disposizione.
In questa applicazione sono state considerate le seguenti funzioni generatrici:
fit.tree_population <- rpart(Sale_Price ~ ., ames_data)
phat.tree_population <- predict(fit.tree_population, newdata = x_test)
fit.rf <- randomForest(Sale_Price ~ ., data = ames_data, ntree = 50, importance = T)
phat.random_forest <- predict(fit.rf, newdata = x_test)
set.seed(123)
ames_bag1 <- bagging(
formula = Sale_Price ~ .,
data = ames_data,
nbagg = 10,
coob = TRUE,
control = rpart.control(minsplit = 2, cp = 0)
)
y_hat_pop_bag<- predict(ames_bag1, newdata = x_test)
Al fine di stimare i modelli sui diversi campioni, si sono implementate le seguenti funzioni che permettono di addestrare un modello e ottenere le previsioni sul test set.
samplePredForDecisionTree <- function(campione){
sample_Tree_Model <- rpart(Sale_Price ~ ., campione)
return(predict(sample_Tree_Model, newdata = x_test))
}
samplePredForRandomForest <- function(campione){
sample_RF_Model <- randomForest(Sale_Price ~ ., data = campione, ntree = 50, importance = T)
return(predict(sample_RF_Model, newdata = x_test))
}
samplePredBagging <- function(campione){
set.seed(123)
library(caret)
library(ipred)
# train bagged model
ames_bag <- bagging(
formula = Sale_Price ~ .,
data = campione,
nbagg = 10,
coob = TRUE,
control = rpart.control(minsplit = 2, cp = 0)
)
return(predict(ames_bag, newdata = x_test))
}
In questa applicazione vengono presi in considerazione 20 modelli per ciascuna ampiezza campionaria. Le ampiezze campionarie di riferimento sono da 100, 300, 500, 700, 900, 1000, 1200 osservazioni.
Con riferimento alla prima ampiezza campionaria formata da 100 osservazioni, si procede costruendo 20 campioni, ciascuno dei quali è ottenuto tramite un campionamento casuale senza reinserimento. Per ogni campione viene addestrato un generico modello, le cui previsioni sono stimate sul test set. Ciò permette di calcolare empiricamente il bias e la varianza del modello.
Si ripete questa procedura per le diverse ampiezze campionarie, in questo modo è possibile valutare l’andamento del bias e varianza all’aumentare della dimensione campionaria.
nModel <-20
noOfSamples = c(100,300,500,700,900,1000,1200)
samplePredictionsTree <- data.frame()
samplePredictionsRF <- data.frame()
samplePredictionsBag <- data.frame()
risultato_finale = list()
risultato_finale_rf = list()
risultato_finale_bg = list()
for (numerosita in noOfSamples) {
for (i in 1:nModel) {
campione = train_data[sample(nrow(ames_train),numerosita, replace = FALSE),]
samplePredictionsTree_temp = samplePredForDecisionTree(campione)
samplePredictionsTree <- rbind(samplePredictionsTree,samplePredictionsTree_temp)
samplePredictionsRF_temp = samplePredForRandomForest(campione)
samplePredictionsRF <- rbind(samplePredictionsRF,samplePredictionsRF_temp)
samplePredictionsBag_temp = samplePredBagging(campione)
samplePredictionsBag <- rbind(samplePredictionsBag,samplePredictionsBag_temp)
}
var_model = calculate_variance_of_model(samplePredictionsTree, y_test)
bias_model = calculate_bias_of_model(samplePredictionsTree, phat.tree_population)
var_model_rf = calculate_variance_of_model(samplePredictionsRF, y_test)
bias_model_rf = calculate_bias_of_model(samplePredictionsRF, phat.random_forest)
var_model_bg = calculate_variance_of_model(samplePredictionsBag, y_test)
bias_model_bg = calculate_bias_of_model(samplePredictionsBag, y_hat_pop_bag)
risultato_finale_bg[[numerosita]] <- c(numerosita,var_model_bg[[1]], var_model_bg[[2]],bias_model_bg)
risultato_finale_rf[[numerosita]] <- c(numerosita,var_model_rf[[1]], var_model_rf[[2]],bias_model_rf)
risultato_finale[[numerosita]] <- c(numerosita,var_model[[1]], var_model[[2]],bias_model)
}
Le seguenti rappresentazioni grafiche mostrano l’andamento di \(bias^2\), \(varianza\) e \(MSE\) al variare dell’ampiezza campionaria. All’aumentare dell’ampiezza campionaria aumenta la \(varianza\), ma si riduce il \(bias^2\), non esiste un modello che minimizza contemporaneamente la distorsione e la varianza, è necessario determinare un “compromesso”.
risultato_finale <- do.call("rbind",risultato_finale)
colnames(risultato_finale) <- c("numerosita","varianza","mse", "bias")
res_tree <- as_tibble(risultato_finale)
res_tree$method <-"Decision Tree"
res_tree
## # A tibble: 7 x 5
## numerosita varianza mse bias method
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 100 0.000705 0.0144 0.00536 Decision Tree
## 2 300 0.00272 0.0136 0.00280 Decision Tree
## 3 500 0.00345 0.0128 0.00217 Decision Tree
## 4 700 0.00344 0.0121 0.00157 Decision Tree
## 5 900 0.00317 0.0118 0.00134 Decision Tree
## 6 1000 0.00308 0.0115 0.00135 Decision Tree
## 7 1200 0.00292 0.0112 0.00131 Decision Tree
require(gridExtra)
plot1 <- ggplot(res_tree, aes(x=numerosita, y=bias))+
geom_line(stat = "identity", col="blue")+
ggtitle("Bias")
plot2 <- ggplot(res_tree, aes(x=numerosita, y=varianza))+
geom_line(stat = "identity", col="red")+
ggtitle("Varianza")
plot3 <-ggplot(res_tree, aes(x=numerosita, y=mse))+
geom_line(stat = "identity", col="green")+
ggtitle("MSE")
grid.arrange(plot1, plot2, plot3, ncol=2)
risultato_finale_rf <- do.call("rbind",risultato_finale_rf)
colnames(risultato_finale_rf) <- c("numerosita","varianza","mse", "bias")
res_rf <- as_tibble(risultato_finale_rf)
res_rf$method <-"Random Forest"
res_rf
## # A tibble: 7 x 5
## numerosita varianza mse bias method
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 100 0.000195 0.00830 0.00257 Random Forest
## 2 300 0.000511 0.00819 0.00223 Random Forest
## 3 500 0.000853 0.00762 0.00184 Random Forest
## 4 700 0.000913 0.00710 0.00157 Random Forest
## 5 900 0.000897 0.00681 0.00144 Random Forest
## 6 1000 0.000913 0.00673 0.00137 Random Forest
## 7 1200 0.000878 0.00654 0.00130 Random Forest
plot1 <- ggplot(res_rf, aes(x=numerosita, y=bias))+
geom_line(stat = "identity", col="blue")+
ggtitle("Bias")
plot2 <- ggplot(res_rf, aes(x=numerosita, y=varianza))+
geom_line(stat = "identity", col="red")+
ggtitle("Varianza")
plot3 <- ggplot(res_rf, aes(x=numerosita, y=mse))+
geom_line(stat = "identity", col="green")+
ggtitle("MSE")
grid.arrange(plot1, plot2, plot3, ncol=2)
risultato_finale_bg <- do.call("rbind",risultato_finale_bg)
colnames(risultato_finale_bg) <- c("numerosita","varianza","mse", "bias")
res_bg <- as_tibble(risultato_finale_bg)
res_bg$method <-"Bagging"
res_bg
## # A tibble: 7 x 5
## numerosita varianza mse bias method
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 100 0.000356 0.00975 0.00354 Bagging
## 2 300 0.00135 0.0101 0.00301 Bagging
## 3 500 0.00195 0.00944 0.00250 Bagging
## 4 700 0.00201 0.00875 0.00213 Bagging
## 5 900 0.00194 0.00845 0.00201 Bagging
## 6 1000 0.00198 0.00833 0.00191 Bagging
## 7 1200 0.00189 0.00803 0.00182 Bagging
plot1 <- ggplot(res_bg, aes(x=numerosita, y=bias))+
geom_line(stat = "identity", col="blue")+
ggtitle("Bias")
plot2 <- ggplot(res_bg, aes(x=numerosita, y=varianza))+
geom_line(stat = "identity", col="red")+
ggtitle("Varianza")
plot3 <- ggplot(res_bg, aes(x=numerosita, y=mse))+
geom_line(stat = "identity", col="green")+
ggtitle("MSE")
grid.arrange(plot1, plot2, plot3, ncol=2)
La seguente rappresentazione confronta \(MSE\), \(bias^2\) e \(varianza\) per i tre modelli considerati fissata un’ampiezza campionaria pari a 1200 osservazioni. Dai risultati emerge che il random forest ha varianza minore, infatti a differenza del bagging, per ciascun nodo anziché esplorare tutti i possibili predittori, viene preso in considerazione solamente un sottoinsieme costituito da \(m<p\) predittori.
res_tot <- rbind(res_tree[7,], res_rf[7,],res_bg[7,])
ggplot(melt(res_tot[,2:5], id.vars='method'),aes(method,value,fill=variable))+
geom_bar(stat="identity",position="dodge")
In generale, in termini di errore quadratico medio il random forest risulta essere il modello migliore.
NB: Il Bagging è un caso particolare di Random Forest, infatti ponendo \(m = p\) si ottengono gli stessi risultati. Nel Random Forest le righe vengono ricampionate nello stesso modo del bagging (ricampionamento con reinserimento) ma viene considerato solamente un sottoinsieme dei predittori. Il difetto principale del Bagging deriva dal limite insito del bootstrap: ovvero, i training set possono essere molto simili tra loro portando così ad ottenere degli alberi molto correlati. Il Random Forest si pone l’obiettivo di rendere gli alberi costruiti dai campioni bootstrap più diversi possibili.
library(caret)
library(glmnet)
library(tidymodels)
library(AmesHousing)
library(hdi)
load("C:/Users/Stefano/Downloads/ames_h_stats (1).RData")
ames <- make_ames()
set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)
set.seed(24873)
ames_folds <- vfold_cv(ames_train)
ames_ind <- rsample2caret(ames_folds)
Si costruisce un modello con variabile target Sale_Price e con 18 variabili esplicative: Bldg_Type, Neighborhood, Year_Built, Gr_Liv_Area, Full_Bath, Year_Sold, Lot_Area, Central_Air, Longitude, Latitude, MS_SubClass, Alley, Lot_Frontage, Pool_Area, Garage_Finish, Foundation, Land_Contour, Roof_Style. La variabile target viene sottoposta a trasformazione logaritmica, mentre le variabili Lot_Area, Gr_Liv_Area, Lot_Frontage vengono normalizzate con una trasformazione di BoxCox. La dimensionalità della variabile Neighborhood viene ridotta includendo in un’unica classe “other” tutte le osservazioni con frequenza relativa <0.05. Tutte le variabili nominali vengono trasformate in variabili dummy e tutte le variabili esplicative vengono normalizzate o eliminate se hanno un solo valore identico in ogni osservazione.
main_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_zv(all_predictors()) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors()) %>%
step_bs(Longitude, Latitude, options = list(df = 5))
main_rec
## Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 18
##
## Operations:
##
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood
## Dummy variables from all_nominal()
## Zero variance filter on all_predictors()
## Centering for all_predictors()
## Scaling for all_predictors()
## B-Splines on Longitude, Latitude
Considerata la distribuzione della statistica H ottenuta a partire dal metodo bootstrap, si sono identificati i potenziali predittori fissando un cutoff pari a 0.001.
int_vars <-
h_stats %>%
dplyr::filter(Estimator == "Bootstrap" & H > 0.001) %>%
pull(Predictor)
Si procede creando un dataframe contenente tutte le possibile coppie di interazioni \(\frac{p(p−1)}{2}\) tra i regressori aventi un valore della statistica H maggiore della soglia.
interactions <- t(combn(as.character(int_vars), 2))
colnames(interactions) <- c("var1", "var2")
interactions
## var1 var2
## [1,] "Garage_Finish" "Roof_Style"
## [2,] "Garage_Finish" "Central_Air"
## [3,] "Garage_Finish" "Year_Built"
## [4,] "Garage_Finish" "Neighborhood"
## [5,] "Garage_Finish" "Gr_Liv_Area"
## [6,] "Roof_Style" "Central_Air"
## [7,] "Roof_Style" "Year_Built"
## [8,] "Roof_Style" "Neighborhood"
## [9,] "Roof_Style" "Gr_Liv_Area"
## [10,] "Central_Air" "Year_Built"
## [11,] "Central_Air" "Neighborhood"
## [12,] "Central_Air" "Gr_Liv_Area"
## [13,] "Year_Built" "Neighborhood"
## [14,] "Year_Built" "Gr_Liv_Area"
## [15,] "Neighborhood" "Gr_Liv_Area"
interactions <-
interactions %>%
as_tibble() %>%
mutate(
term =
paste0(
"starts_with('",
var1,
"'):starts_with('",
var2,
"')"
)
) %>%
pull(term) %>%
paste(collapse = "+")
interactions <- paste("~", interactions)
interactions <- as.formula(interactions)
int_rec <-
recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
Central_Air + Longitude + Latitude + MS_SubClass +
Alley + Lot_Frontage + Pool_Area + Garage_Finish +
Foundation + Land_Contour + Roof_Style,
data = ames_train) %>%
step_log(Sale_Price, base = 10) %>%
step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
step_other(Neighborhood, threshold = 0.05) %>%
step_dummy(all_nominal()) %>%
step_interact(interactions) %>%
step_zv(all_predictors()) %>%
step_center(all_predictors()) %>%
step_scale(all_predictors()) %>%
step_bs(Longitude, Latitude, options = list(df = 5))
La ricerca dei valori ottimi di \(\alpha\) e \(\lambda\) avviene in una griglia contenente i possibili valori per tali parametri.
ctrl <-
trainControl(
method = "cv",
index = ames_ind$index,
indexOut = ames_ind$indexOut
)
glmn_grid <- expand.grid(alpha = seq(.2, 1, by = .2), lambda = 10^seq(-4, -1, by = 0.1))
main_glmn_h <-
train(main_rec,
data = ames_train,
method = "glmnet",
tuneGrid = glmn_grid,
trControl = ctrl
)
int_glmn_h <-
train(int_rec,
data = ames_train,
method = "glmnet",
tuneGrid = glmn_grid,
trControl = ctrl
)
int_glmn_h$bestTune
## alpha lambda
## 136 1 0.001258925
p <-
ggplot(int_glmn_h) +
scale_x_log10() +
theme_bw() +
theme(legend.position = "top")
p
Questo grafico mostra l’andamento dell’RMSE al variare dei parametri \(\lambda\) e \(\alpha\). Ogni curva corrisponde ad un diverso valore di \(\alpha\) e rappresenta le performance del modello stimato in funzione del parametro di penalizzazione \(\lambda\). In questo modo è possibile individuare i valori dei parametri di tuning che minimizzano l’RMSE.
Il nostro obiettivo è quello di selezionare tra i 157 regressori quelli che risultano maggiormente implicati nello spiegare la variazione della variabile dipendente. La ricerca di lambda ottima avviene nella griglia di valori specificata all’interno del comando cv.glmnet.
lasso = cv.glmnet(as.matrix(train_data[,-8]),
train_data$Sale_Price, alpha=1,
lambda = c(10^seq(-4, -1, by = 0.1)))
plot(lasso, ylim = c(0,0.04))
In questo grafico si evince il lambda path: andamento della curva dell’errore quadratico medio di previsione dato dalla cross validation. Per ciascun valore di lambda la media dell’errore di cross validation è rappresentata dal punto rosso le barre rappresentano gli standard error della media in corrispodenza di quel valore di lambda.
##
## Call: cv.glmnet(x = as.matrix(train_data[, -8]), y = train_data$Sale_Price, lambda = c(10^seq(-4, -1, by = 0.1)), alpha = 1)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.001259 20 0.005689 0.0005258 67
## 1se 0.005012 14 0.006049 0.0005793 32
Il valore di \(\lambda\) che minimizza la curva di cross validation permette di stabilire che il numero di predittori con coefficente di regressione diverso da 0 è 67.
NB.
É possibile “tollerare” un piccolo incremento nell errore guadagnandoci dal punto di vista interpretativo; in questo modo scelgo un nuovo valore di \(\lambda_{1se}\) che permette di ottenere un errore medio di previsione maggiore ma al contempo restituisce un modello più interpretabile.
Per valutare l’effetto delle interazioni sulla capacità di prevedere la varibile dipendente, si confronta l’errore quadratico medio dei due modelli sul test set.
yhats_main = predict(main_glmn_h, ames_test)
yhats_int = predict(int_glmn_h, ames_test)
## [1] "Modello main: 0.0053"
## [1] "Modello main + interazioni: 0.0054"
In termini di errore quadratico medio i due modelli risultano pressoché uguali, nonostante questo si predilige il modello senza interazioni perchè meno complesso e più facilmente interpretabile.
Considerando che la funzione generatrice dei dati è ignota, si assume che questa coincida con il modello stimato sulla totalità dei dati a nostra disposizione.
In questa applicazione è stata considerata la seguente funzione generatrice, che coincide con il modello lasso caratterizzato dal valore di \(\lambda\) che minimizza l’errore quadratico medio.
pop_lasso_model <- glmnet(ames_data[,-8], ames_data$Sale_Price, alpha = 1, lambda = 0.001258925)
y_hat_pop_lasso <- predict(pop_lasso_model, newx = as.matrix(x_test))
La seguente funzione consente di addestrare il modello e ottenere le previsioni sul test set.
samplePredForLasso <- function(campione, lambda){
sample_Lasso_Model <- glmnet(campione[,-8], campione$Sale_Price, alpha = 1, lambda = lambda)
return(predict(sample_Lasso_Model, newx = as.matrix(x_test)))
}
In questa applicazione vengono stimati 30 modelli con un’ampiezza campionaria pari a 2000, ciascuno dei quali è ottenuto tramite un campionamento casuale con reinserimento.
Fissato un valore di \(\lambda\), per ogni campione viene stimato un modello lasso e si effettuano previsioni sul test set. Ciò permette di calcolare empiricamente il bias e la varianza del modello.
Si ripete questa procedura per i diversi valori di \(\lambda\), in questo modo è possibile valutare l’andamento del bias e varianza all’aumentare della dimensione campionaria.
nModel <-30
samplePredictionsLasso <- data.frame()
risultato_finale = list()
numerosita=2000
lambdas = c(0.001258925, 0.01, 0.05,0.1, 0.15, 0.2, 0.25, 0.30, 0.35, 0.40)
j=1
for (lambda in lambdas) {
for (i in 1:nModel) {
campione = train_data[sample(nrow(ames_train),numerosita, replace = TRUE),]
samplePredictionsLasso_temp = samplePredForLasso(campione, lambda)
samplePredictionsLasso <- rbind(samplePredictionsLasso,samplePredictionsLasso_temp)
}
var_model = calculate_varaince_of_model(samplePredictionsLasso, y_test)#aggiungere mse
bias_model = calculate_bias_of_model(samplePredictionsLasso, y_hat_pop_lasso)#aggiungere mse
risultato_finale[[j]] <- c(lambda,var_model[[1]], var_model[[2]],bias_model)
j=j+1
}
risultato_finale <- do.call("rbind",risultato_finale)
colnames(risultato_finale) <- c("lambda","varianza","mse", "bias")
risultato_finale <- as_tibble(risultato_finale)
risultato_finale
## # A tibble: 10 x 4
## lambda varianza mse bias
## <dbl> <dbl> <dbl> <dbl>
## 1 0.00126 0.0254 0.0386 0.0168
## 2 0.01 0.0229 0.0363 0.0168
## 3 0.05 0.0183 0.0316 0.0168
## 4 0.1 0.0142 0.0273 0.0168
## 5 0.15 0.0114 0.0244 0.0168
## 6 0.2 0.00949 0.0224 0.0168
## 7 0.25 0.00814 0.0209 0.0168
## 8 0.3 0.00712 0.0198 0.0168
## 9 0.35 0.00634 0.0190 0.0169
## 10 0.4 0.00570 0.0183 0.0169
ggplot(risultato_finale, aes(x=lambda, y=bias))+
geom_line(stat = "identity", col="blue")+
ggtitle("Bias")
ggplot(risultato_finale, aes(x=lambda, y=varianza))+
geom_line(stat = "identity", col="red")+
ggtitle("Varianza")
ggplot(risultato_finale, aes(x=lambda, y=mse))+
geom_line(stat = "identity", col="green")+
ggtitle("MSE")
Da questi grafici si evince la presenza di un trade-off tra \(distorsione^2\) e varianza, in quanto all’aumentare della penalizzazione diminuisce la varianza, ma aumenta il \(bias^2\); questo perchè il modello lasso esclude un numero sempre maggiore di predittori.